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ABSTRACT 

Following Bi & Davidsen (1997), we perform one dimensional semi analytic 
simulations along the lines of sight to model the intergalactic medium (IGM). 
Since this procedure is computationally efficient in probing the parameter space 
- and reasonably accurate - we use it to recover the values of various parameters 
related to the IGM (for a fixed background cosmology) by comparing the model 
predictions with different observations. For the currently favoured LCDM model 
(Q m = 0.4, Q\ = 0.6 and h = 0.65), we obtain, using statistics obtained from the 
transmitted flux, constraints on (i) the combination / = (£7#/i 2 ) 2 / J-12, where 
Qb is the baryonic density parameter and J-12 is the total photoionisation rate 
in units of 10~ 12 s -1 , (ii) temperature T corresponding to the mean density and 
(iii) the slope 7 of the effective equation of state of the IGM at a mean redshift 
z ~ 2.5. We find that 0.8 < (T /10 4 K) < 2.5 and 1.3 < 7 < 2.3. while the 
constraint obtained on / is 0.020 2 < f < 0.032 2 . A reliable lower bound on J_ 12 
can be used to put a lower bound on f2#/i 2 , which can be compared with similar 
constraints obtained from Big Bang Nucleosynthesis (BBN) and CMBR studies. 
We find that if J_i2 > 1.2, the lower bound on Qsh 2 is in violation of the BBN 
value. 

Subject headings: cosmology: large-scale structure of universe - intergalactic 
medium - quasars: absorption lines 
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1. Introduction 

A significant fraction of the baryons at z < 5 are found in the form of a diffuse inter- 
galactic medium (IGM), which is usually probed through the absorption lines produced by 
them on the spectrum of the distant QSOs. It is believed that while the metal line systems 
(detected through Mg II or C IV doublets) seen in the QSO spectra could be associated 
with the halos of the intervening luminous galaxies (Bergeron & Boisse 1991; Steidel 1993), 
most of the low neutral hydrogen column density absorption lines (commonly called as 'Lya' 
clouds) are due to the low amplitude baryonic fluctuations in the IGM. 

Probing the baryonic structure formation through Lya absorption lines has two advan- 
tages. First, there are large number of absorption lines. Typically, one can observe more 
than a few hundreds of lines per unit redshift range along any one line of sight. This pro- 
vides us with a large unbiased dataset, using which the statistical studies can be performed 
efficiently The second advantage is that the Lya absorption lines are more straightforward 
to model than, say, luminous galaxies. The modelling of galaxies is complicated by the fact 
that one has to take into account processes like the star formation, radiation feedback and 
so on; these processes are not that effective in the IGM, and one can ignore them at the first 
approximation. 

The study of the IGM can - potentially - provide us with information about different 
aspects of the the baryonic structures the universe like (i) the mass power spectrum (Croft 
et al. 1998; Hui 1999; Croft et al. 1999), (ii) the total baryonic density (Sl B ) and the 
total photoionisation rate due to the local ionising background radiation (J) and (iii) the 
reionisation history of the universe (Hui & Gnedin 1997). 

There have been various numerical and semi analytical models in the literature for the 
IGM, all of which are based on the view that the Lya clouds are small scale density fluctu- 
ations as predicted by the models of structure formation. The hydrodynamical simulations 
(Bond, Szalay & Silk 1988; Cen et al. 1994; Zhang, Anninos & Norman 1995; Hernquist et 
al. 1996; Miralda-Escude et al. 1996; Riediger, Petitjean & Miicket 1998; Theuns, Leonard 
& Efstathiou 1998; Theuns et al. 1998; Dave et al. 1999) incorporate most of the ongoing 
physical processes in the IGM and hence they are necessary for understanding the evolution 
of the IGM. However, due to limited numerical resolution and computing power, they are 
able to probe only a small box size (10-20 Mpc). Hence people have tried to complement 
the numerical studies with analytic and semi analytic ones (Doroshkevich & Shandarin 1977; 
McGill 1990; Bi 1993; Bi, Ge & Fang 1995; Gnedin & Hui 1996; Hui, Gnedin & Zhang 1997; 
Bi & Davidsen 1997, hereafter BD; Choudhury, Padmanabhan & Srianand 2000, hereafter 
Paper I). The semi analytic models do not have problems related to limited numerical reso- 
lution or box sizes, and can be used to probe a wide range of parameters. 
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The numerical simulations suggest that most of the Lya lines arise due to linear or 
quasi-linear density fluctuations. Therefore one can neglect the highly non-linear baryonic 
processes (like shock heating) as a first approximation. However, since a simple linear density 
evolution cannot produce the saturated Lya systems, one cannot completely ignore the non- 
linear effects. The non-linear baryonic density can be calculated from the linear one using 
some approximation scheme, like the Zeldovich approximation (Doroshkevich & Shandarin 
1977; McGill 1990; Hui, Gnedin & Zhang 1997), or the lognormal approximation (Bi 1993; 
Gnedin & Hui 1996; BD; Paper I). The neutral fraction is then estimated by considering the 
equilibrium between the rate of photoionisation due to background radiation and the rate 
of recombination estimated from the temperature defined through the equation of state. All 
these models depend on various IGM parameters such as Qb, J, equation of state and the 
Jeans length, as well as the cosmological parameters like Q m , etc. 

Besides using these simplifying assumptions, BD realised that it is sufficient to simulate 
the IGM in ID rather than in 3D. This increases the computing power drastically, and one 
can probe large box sizes (hundreds of Mpc) with high enough resolution. BD performed a 
detailed study of the evolution of the IGM from z = 2 to 4. They also compared their pre- 
dictions of column density distribution with hydrodynamical simulations and observations, 
and found a good agreement. In this work, we follow the idea proposed by BD and carry out 
semi analytic simulations of the low density IGM. The results obtained from such simulations 
are found to be in quite good agreement with various observations (as described in Section 
5), thus indicating that the lognormal approximation might be a reasonable assumption for 
the low density IGM. We extend our studies to probe the parameter space and constrain 
the parameters for a particular redshift bin using different statistics obtained from the spec- 
trum. Since the recovery of cosmological parameters is not possible with ill constrained IGM 
parameters (for a detailed discussion, see Paper I), we concentrate only on the parameters 
related to the IGM at a particular redshift (in this case, z = 2.41). The parameters are 
the slope of the equation of state (7), the temperature corresponding to the mean baryonic 
density or the mean temperature(T ) [in this paper, we shall use the term 'mean tempera- 
ture' to be equivalent to the temperature corresponding to the mean baryonic density] and 
a combination of the baryonic density parameter (Qb) an d the total photoionisation rate 
due to the local ionising radiation field (J), the combination being / = (Qsh 2 ) 2 / J-12, where 
J_ 12 = .//(lO^V 1 ). We find that different statist ics are sensitive to different parameters, 
and hence they can be used simultaneously to constrain the parameter space. 

In previous studies involving numerical simulations, the parameter / is usually deter- 
mined by demanding that the simulated mean transmitted flux match with the observations 
(Rauch et al. 1997; McDonald et al. 2000a). Then for a given /, the constraints on 7 and 
T are usually obtained by fitting the lower envelope of the Nm — b scatter plot (Schaye 



-4- 



et al. 2000; Ricotti, Gnedin & Shull 2000; Bryan k Machacek 2000; McDonald et al. 
2000b), where Nm is the column density and b is the thermal velocity dispersion (defined as 
b = ^Ik^T jm p ). However, because of limited box size in the hydrodynamical simulations, 
the continuum of the transmitted flux is not well identified and this introduces errors in 
the calculation of the mean transmitted flux. In our approach, we constrain all the three 
parameters simultaneously using all the available transmitted flux statistics, thus utilising 
all the information available in the spectrum. 

Section 2 gives the basic structure of the simulation strategy. Although the basic idea is 
the same as in BD, we repeat some of the details for completeness. This section also discusses 
about the various assumptions used at different stages. We discuss the various parameters 
used to model the simulation in Section 3. Section 4 contains a very brief discussion on the 
various statistical quantities studied in this paper. The next section contains the results, 
where we compare our simulations with available observational data and constrain /, 7 and 
T . Finally, we summarise our conclusions in section 6. 



2. Basic Outline of the Simulation 

We describe the structure of the numerical simulation in this section, which is essentially 
the same as in BD, for completeness and setting up the notation. 

Let -Pdm(^) denote the linear DM power spectrum in 3D at the present epoch (z — 0). 
Then the power spectrum for any arbitrary z is given by 

p{V(k,z) = D*(z)P^(k), (1) 

where D(z) gives the evolution of the linear density contrast. The linear baryonic power 
spectrum is related to the DM power spectrum through the relation (Fang et al. 1993) 

p(3)( k z ) _ P Dm(M) (2] 



where 

x b(z) = — 

-"0 



2<yk B T„ 



1/2 

(3) 

3fim p il m (l + z) J 

is the Jeans length; ji is the mean molecular weight of the IGM, given by /z = 4/(8 — 5Y), 
where Y is the helium weight fraction. (This relation assumes that the IGM consists mostly 
of fully ionised hydrogen and helium. In this paper, we take Y = 0.24.) T rn is the density 
averaged temperature of the IGM and 7 is the ratio of specific heats. il m is the cosmological 
density parameter. Strictly speaking, equation (2) is valid only for the case where Xf, is 
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independent of z, but it is shown by Bi, Borner & Chu (1992) that equation (2) is a good 
approximation for P B 3 \k, z) even when Xb has a redshift dependence. 

At this point, it is appropriate to stress some features of the parameter T m . The obvious 
interpretation of T m will be as the mean temperature of the IGM, T (the temperature at 
the mean density), i.e., T m = T . However, according to BD, using T in equation (3) 
leads to a value of the linear baryonic density fluctuation, <tb, larger than what we expect 
from hydrodynamical simulations. Hence, they suggested the use of a density averaged 
temperature. Since T m appears only in the expression for the Jeans length (equation (3)), it 
can also be defined as the effective temperature which determines the Jeans length. It is clear 
that the combination 7T m can, in principle, be fixed if erg is known through hydrodynamical 
simulations. In this work, we choose a B (z = 2.41) to be 1.34, which gives 7T m = 5.115 x 
10 4 K. Our choice of ub(z = 2.41) is consistent with that of BD (see their Figure 3). Also, 
simulations of Carlberg & Couchman (1989) give ub(z = 2.8) = 0.953, but one should note 
that the power spectrum they used was normalised to a value that was 1.4 times smaller than 
ours. Furthermore, using a power spectrum normalised to a value 1.3 times larger than ours, 
Gnedin (1998) obtains &b{z = 2.85) = 2.25. All these values are consistent with our choice 
for which the value of Jeans length at z — 2.41 is O.^Q^^h -1 Mpc. For the background 
cosmology with Vt m = 0.4, = 0.6, this corresponds to a velocity scale of 22.3 km s _1 . 
Since we shall be mainly concerned with a small redshift bin (Az = 0.58), the evolution of 
T m should not affect the results significantly. Hence, we take T m to be independent of z. 

Once the power spectrum of linear density perturbations in 3D is obtained, one can 
obtain the corresponding power spectra for density (as well as velocity) perturbations in ID. 
One can show that the baryonic power spectrum in ID is given by 




1*1 



(4) 



while the power spectrum for linear velocity perturbations in ID is 





where a is the scale factor and a is given by the Friedman equations 





with 




(7) 
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The density and the velocity fields are correlated with the correlation being given by 

1 f°° Ah 1 

Pj£{k,z) = i ^r p B\k',z). (8) 

To simulate the density and velocity fields in ID for a particular redshift z, we follow 
the procedure given by Bi (1993). We start with two independent Gaussian fields, w (k) and 
uo(k), having unit power spectrum, i.e., (wQ(k)wo(p)) = (uQ(k)uo(p)) = 27i5j)irac(k — p). We 
can then get two independent Gaussian fields having power spectra P w (k,z) and P u (k,z) 
respectively 

w(k, z) — w (k) yj P w (k, z) , u(k, z)=U (k)y/P u (k,z). (9) 
We choose these power spectra to be of the following form (Bi 1993; BD) 

1 /-OO JK/ 

P w (k,z)=(3- 1 (k,z)- -rrP ( B\k',z) (10) 
2tt J ]k] k' 

and 

P u (k,z) = — dk' k' P B 3 \k',z)-P w (k,z), (11) 

^ J\k\ 

where 

(3(k, z) = J{k{K ' } ^ \ (12) 
f~(dk'/k>)Pg\k',z) 

The linear density and the velocity fields in the /c-space are then given by 

Ssik, z) — w(k, z) + u(k, z), (13) 

v(k, z) = i ak[3(k, z)w(k, z). (14) 

The corresponding fields Sb(x,z) and v(x,z) in the real comoving space are obtained by 
using Fourier transforms. One should keep in mind that the above analysis is done in the 
framework of linear perturbation theory. 

However, to study the properties of the IGM one has to take into account the non- 
linearities in the density distribution and various physical processes such as shocks, radiation 
field, cooling etc. Detailed hydrodynamical modelling of IGM has shown that most of the 
low column density Lya absorption (i.e. iVm < 10 14 cm~ 2 ) are produced by regions that 
are either in the linear or in the weakly non-linear regime (Cen et al. 1994; Zhang et al. 
1995; Hernquist et al. 1996; Miralda-Escude et al. 1996; Theuns, Leonard & Efstathiou 
1998; Theuns et al. 1998; Dave et al. 1999). The lower envelope of the column density, 
Nm — b scatter plot (Schaye et al. 1999; Schaye et al. 2000) suggests that there is a well 
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defined relationship between the density and the temperature of the IGM (Hui & Gnedin 
1997). Thus it is possible to model low column density systems using simple prescription for 
the non-linear density field and an equation of state. 

Following BD, we take into account the effect of non-linearities of density perturbations 
by assuming the number density distribution of the baryons, ub(x,z) to be a lognormal 
random field 

n B (x, z) = A e SB( - x ^ (15) 

where 5b(x, z) is the linear density contrast in baryons, and A is a constant to be determined. 
The mean value of ub(x, z) is given by 

{n B (x,z))=n (z)=A(e 5B ^), (16) 

where no(z) is related to the baryonic density parameter Q B through the relation 

n Q {z) = ^{l + zf. (17) 
fi B m p 

Here p c = 1.8791 x I0~ 29 h 2 cm~ 3 is the critical density of the universe and [i B m v is the 
mass per baryonic particle, given by /xb7b p = Am p / (4 — 3Y). Hence, we get the value of the 
constant as 

U§ B (x,z)\ v i 

and 

n B (x,z)=n (z) {eSB(x<z) y (19) 

The lognormal distribution was introduced by Coles & Jones (1991) as a model for the non- 
linear matter distribution in the universe. Detailed arguments as to why this ansatz should 
be reasonable in studying non-linear density distribution can be found in Coles & Jones 
(1991), BD and Paper I. However, we would still like to stress some points regarding the use 
of the lognormal ansatz for baryons in the current context. 

In the past, there have been attempts to use the lognormal distribution to model the dark 
matter. However, we now know (based on Non-linear Scaling Relations; see Nityananda & 
Padmanabhan 1994; Padmanabhan 1996) that any local mapping of the form 5nl = F[8l] 
is bad for dark matter (also see Coles, Melott & Shandarin 1993). There is, however, a 
strong theoretical argument (see Paper I) which shows that lognormal produces the correct 
limits at the two extremes for baryons. At large spatial scales, where the density contrast 
is small (5b < 1), equation (19) reduces to ne/no — 1 + 5b, which is just what we expect 
from linear theory. More importantly, on small scales, equation (19) becomes the isothermal 
hydrostatic solution, which describes highly clumped structures like intracluster gas, ne oc 
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exp(—fj,m p i/jBM/jkBT), where -0dm is the dark matter potential (Sarazin & Bahcall 1977). 
This gives an indication that even though the lognormal ansatz is poor for dark matter 
distribution, it might still work reasonably well for baryonic matter since it is correctly 
constrained at both extremes. 

Comparison with full hydrodynamical simulations reinforces this conclusion. BD have 
used the results from the hydrodynamical simulations of Miralda-Escude et al. (1996), and 
found that the baryonic density distribution can be well fitted by a lognormal function at 
z = 3. Also, the range of parameters for the IGM recovered by us and those by the full 
hydrodynamical simulations of McDonald et al. (2000a, b) agree quite well (as discussed later 
in Section 6). This shows that the lognormal assumption agrees with the hydrodynamical 
simulations as well. 

The above arguments should convince the reader that the lognormal assumption, in spite 
of its limitations, provides us with a tool in studying the baryonic structure formation semi 
analytically. No such approximation can reproduce the results obtained from the full hy- 
drodynamical simulations exactly and their values lie in providing faster route to reasonably 
accurate results. Our attempts should be viewed in the backdrop of such a philosophy. 

There is no obvious way to deal with the non-linearities in the velocity field, but fortu- 
nately this is not needed in the current work as can be seen from the following argument. 
The velocity field plays two separate roles in the context of our work. The first one is that 
the velocity determines the movement of the individual particles at a given instant of time, 
which in turn affects the underlying density field in the next time step. Mathematically this 
feature is represented by the Euler equation, which connects the density field to velocity field. 
Given any prescription for density field, the Euler equation implicitly leads to a consistent 
velocity field. Hence this dynamical effect of velocity field - viz. moving the mass to the 
right location - is indirectly taken into account in any prescription for non-linear density. In 
our case, the lognormal ansatz takes care of this feature. The second effect of the velocity 
field is purely kinematic - it shifts the positions of the absorption lines. In our work we 
will be using the thermal velocity dispersion b in the Voigt profile while analysing the lines. 
Since the Lja absorption lines originate from quasi-linear density regions, the velocity field 
will be subdominant or of the same order as the thermal velocity dispersion b, which will be 
taken care of in our analysis. 

Once the non-linear baryon density is obtained, it is trivial to get the fraction of hydrogen 
in the neutral form, / HI , in the IGM by solving the ionisation equilibrium equation for 
hydrogen 

a(T)n p n e = T ci (T)n e n m + Jn m , (20) 
where a(T) is the radiative recombination rate, r ci (T) is the rate of collisional ionisation and 
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J is rate of photoionisation for hydrogen (Black 1981); n p , n e and rim are the number densities 
of proton, electron and neutral hydrogen, respectively. In general, all these quantities are 
functions of z and all except J depend on the position x too. We shall parametrise J(z) by 
a dimensionless quantity J_i2(z), defined by J(z) = J_ 12 (-2)10" 12 s _1 . For comparison, we 
mention that our J is equal to the quantity J21G1 used by BD. 



(21) 



Black (1981) gives the approximate form of the recombination and ionisation rates as 
follows: 

q(T) f 4.36 x io- 10 T-°- 7573 (if T > 5000K) 
cm^- 1 ~ \ 2.17 x io- 10 T-°- 6756 (if T < 5000K) 

and 

T ci (T) = 5.85 x 10 _11 T 1/2 exp(-157809.1/T) cmV 1 , (22) 

where T is in Kelvin. One can see that the expression for a(T) diverges as T — > which 
needs to be regularised by a temperature cutoff at the lower end in numerical work. BD have 
used the photoionisation temperature as the minimum temperature, which is about 10 4 K 
(see also Theuns et al. 1998). In situations where T > 10 4 K, we too shall use the same 
value. However, when T < 10 4 K, we have taken the minimum temperature to be 5000K. 

The IGM contains mainly hydrogen, a smaller amount of helium (weight fraction, Y ~ 
0.24) and negligible amount of other heavier elements. In that case we can write 
where n is a constant, greater than but very close to unity. However, in the following 
calculation we have neglected the presence of the heavier elements completely for simplicity. 
Let us define the neutral fraction of hydrogen, / HI by 

rtB n m + n p 

(we ignore the number density contributed be helium because, usually, nne/^B < 0.1). Using 
equation (20) in (23), one gets 

( \ a(T(x,z)) 

f m { X , z) = z)) - Vci{T ^ z)) - J{z)/rie{Xi z) ■ y ) 

We express n e in terms of % by assuming that fm <^ 1 and all the helium present is in the 
fully ionised form. In such case, 

n e /n B = fi e = 2(2 -Y)/ (A -3Y). (25) 

Then, 

/ n a(T(x,z))n B (x,z) 

HU ' ] a(T(x : z)) + r ci (T(x,z)) + J(z)/(fi e n B (x,z)) 1 ) 
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We calculate T(x, z) using the polytropic equation of state 

'n B {x,z) x ' 



T(x,z) = T (z) 



n (z) 



(27) 



where T (z) is the temperature of the IGM at the mean density. We now know the neutral 
hydrogen density at any redshift z along a particular axis. Our next goal is to find the 
density along a line of sight. This can be done by obtaining the density field along the 
backward light cone. In other words, we must obtain the quantity nm(x, z(x)), where x and 
z are related through the expression 



x(z) 



cIh(z') dz' , 



(28) 



with 



d H (z) 



= —[n A + n m (i + zf + n k (i + z) 2 ]-^ 2 . 

-no 



(29) 



Similarly one can also get the velocity field v{x, z(x)) along the same LOS. Once the neutral 
hydrogen density and the velocity along the LOS is known, the Lya absorption optical depth 
at redshift zq can be obtained from the relation 



T ( z o) 



c/. 



1/ 



dx 



n m (x,z(x)) 



b(x,z(x))(l + z(x)) 



x V 



a, 



c(z{x) - z ) 



v(x, z(x)) 

b(x, z(x))(l + zq) ' b(x,z(x)) 



+ 



where 



b(x, z(x)) 



t2k B T(x,z(x)) 



(30) 



(31) 



4.45 x 10 18 cnr and V is the Voigt function. I a is related to the Lya absorption cross 



section through 



<T a (v) 



Clg 

b^ 



V 



a, 



c [v - v n 



bv n 



(32) 



For low column density regions, the natural broadening is not that important, and the Voigt 
function reduces to a simple Gaussian 



r A« / 

M a >~p - ex P ( " 



b 2 



(33) 



Since we are mostly dealing with weakly non-linear regimes, where the densities are not too 
high, this approximation does not introduce any significant error in the final results. The 
optical depths obtained above are used to get the final line of sight spectrum. 
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3. 



Model Parameters 



Following Paper I, the model parameters can be broadly divided into two classes, namely, 
those related to the background cosmology and those related to the baryonic IGM. In Paper 
I, we have considered various CDM cosmological models (SCDM, LCDM, OCDM) and a 
range of IGM parameters, and have found that the parameters related to the background 
cosmology cannot be constrained uniquely with ill defined IGM parameters. Consequently, 
the approach taken in this paper is to use the available observations for constraining the IGM 
related parameters under the framework of most favoured structure formation scenario. We 
consider the following parameters for the LCDM model. The CDM power spectrum in 3D 
is taken to be (Efstathiou, Bond & White 1992) 



where v = 1.13, a x = (6.4/r)^ 1 Mpc, a 2 = (3.0/r)^ 1 Mpc, a 3 = {1.7 /T)^ 1 Mpc and 
T = Vt rn h. The normalisation parameter Adm is fixed through the value of as (the rms 
density fluctuation in spheres of radius 8 hr l Mpc) which is taken to be a 8 = 0.79. The 
other model parameters are: 



The values of Q m and are consistent with the best fitted parameters of Ostriker & 
Steinhardt (1995). The value of <j 8 is obtained from the first year COBE normalisation 
(Kofman, Gnedin & Bahcall 1993). This value is also consistent with those obtained from 
the observed local abundance of clusters by Eke, Cole & Frenk (1996). An identical LCDM 
model is considered in the hydrodynamical simulations by Miralda-Escude et al. (1996). 

Once the cosmology is fixed, we turn our attention towards the parameters related to 
the baryons. 

1. Slope of the effective equation of state (7): It is known that the value of 7, at any 
given epoch, depends on the reionisation history of the universe (Theuns et al. 1998, 
Hui & Gnedin 1997). The value of 7 and its evolution are still quite uncertain. Using 
Voigt profile fits to the observed Lja absorption lines one can in principle obtain the 
value of 7. In this work, we will keep 7 as a free parameter and ignore its redshift 
evolution. 

2. Mean temperature (T (z)): Mean temperature of the IGM is decided by the various 
heating and cooling processes. In addition to the reionisation history local radiation 
field will also affect the value of T (z). In the case of full hydrodynamical models, 




V 



(34) 



tt m = 0.4, Q\ = 0.6, h = 0.65. 



(35) 
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the mean temperature is estimated self-consistently by considering various processes. 
However in our approach we consider the mean temperature as a free parameter. We 
also take it as independent of z within the small redshift bin we consider. 

3. Qsh 2 and J-uiz): If one compares the typical values of the three quantities in the 
denominator in the right hand side of equation (26), one can verify that ct(T) and r ci (T) 
are much smaller compared to J/ub (typically, for T ~ 5 x 10 4 K, a ~ 10 _13 cm 3 s _1 , 
r ci ~ 10~ 14 cm 3 s _1 and J(z)/ub ~ 10~ 5 cm 3 s _1 ). This means that we can write nm — 
an 2 B /J. Since ub oc tiq oc Vtsh 2 and J(z) oc J-u(z), we see that only the combination 
f(z) = (VLsli 2 ) 2 j J-u{z) appears in the expression for optical depth. We shall treat 
this quantity f(z) as a free parameter. We shall also assume that the photoionisation 
rate does not depend on z (at the least, it does not vary considerably within the small 
redshift bin we are interested in). This will make / independent of z. 

So, we finally end up with three free parameters, namely 7, T , and /. 

4. Statistical Quantities: Definitions 

We perform various statistics on our simulated spectrum, as one usually does with the 
real data, to constrain various parameters of our model. 

^From the spectrum, one can immediately calculate the mean transmitted flux (F) and 
the rms flux fluctuations (a 2 ?)- The transmitted flux data can also be used to obtain three 
important statistics (McDonald et al. 2000a). These are: (i) the probability distribution 
function (PDF) for the transmitted flux, (ii) the correlation function of the transmitted 
flux, defined as £(At>) = ((F(v) — F)(F(v + Av) — F)) and (iii) the flux power spectrum 
(Pp(k)). The power spectrum is calculated using the Lomb periodogram technique (Lomb 
1976; Scargle 1982; Press et al. 1992) and the normalisation used is the same as mentioned 
in McDonald et al. (2000a), i.e., a 2 F = ff° (dk/2ir)PF(k). The advantage in the case of 
statistics obtained from transmitted flux is that the numerical procedure is quite fast. We 
use these statistics for constraining the parameter space. 

For a set of most favourable values of parameters, we decompose the spectrum into 
individual lines using Voigt profile analysis, and use them to check our predictions with the 
observations. The statistics used for this purpose are (i) the number of lines (absorbers) per 
unit redshift range (dN/dz) and the column density distribution (f(N m )), defined as the 
number of lines (absorbers) per unit redshift path per unit column density range (Kim et 
al. 1997), (ii) the distribution of the b parameter and (iii) the two point correlation function 
for the absorbers, defined as £ c i oud (Ai;) = [iV sim (At> )/N cxp (Av )] — 1, where N sim (Av) is the 
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number of cloud pairs with a velocity separation Av obtained from the simulated data and 
iV b s (Af ) is the number of pairs expected from a random distribution of clouds (Sargent et 
al. 1980; Webb 1987; Srianand & Khare 1994; Kulkarni et al. 1996; Srianand 1996; Khare 
et al. 1997; Kim et al. 1997; Cristiani et al. 1997). 

5. Results 

In this paper, we have concentrated our studies in the redshift range 2.09-2.67. This 
range corresponds to a box size of 436 h~ x Mpc (for the cosmology we are considering), which 
is virtually impossible to probe in a full 3D hydrodynamical simulation with high enough 
resolution. The number of grid points used in this work along the line of sight was 2 15 , which 
were equispaced in the comoving coordinate x. The simulated flux data was then resampled 
with AA=0.04 A and a random noise of (S/N)=30 was added, exactly as is done with the 
observed data. We mention here that even if we increase the number of points (i.e., try 
to achieve a better resolution), the resampling mentioned above would make sure that the 
statistics obtained from the transmitted flux are not affected. We found that the continuum 
of the spectrum is quite well defined at this redshift, and hence it was not necessary to make 
any extra normalisation. 

5.1. Comparison with observations 

We use various statistics obtained from the observational data given by McDonald et al 
(2000a). In the redshift range 2.09-2.67, they have considered data from 5 QSOs, namely, 
Q2343+123 {z em = 2.52), Q1442+293 (z em = 2.67), KP77: 1623+2653 (z cm = 2.526), 
Q 1107+485 (z em = 3.00) and Q1425+604 (z cm = 3.20). Each of these quasar sight lines 
span different regions of the redshift interval, and hence all the redshift s are not equally 
weighted in the above mentioned redshift range. However in our simulated data we cover 
the same redshift range giving equal weightage. We have confirmed that the correction intro- 
duced due to this uneven weightage in observed data is negligible (i.e., much below typical 
observational errors). 

The allowed range for various parameters are obtained by demanding that the simulated 
data pass through most of the observed points, within the allowed la error limits. The value 
of / is strongly constrained between 0.020 2 -0.032 2 (regardless of the value of 7 and T ). 
In the above range of /, we consider models with / = 0.023 2 , 0.026 2 , 0.029 2 and obtain 
constrains on T and 7 so that all the observed statistics obtained from the transmitted flux 
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f=0.029 3 f = 0.026 3 f = 0.023 2 




1.5 2 2.5 1.5 2 2.5 1.5 2 2.5 
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Fig. 1. — The constraints obtained in the 7 — To space for different values of /, using 
transmitted flux statistics. The shaded regions denote the range allowed by observations. 
The boundaries are uncertain by an amount 0.1 along 7 axis and by 1000K along the T axis 
because of finite sampling, which is shown by a cross at the upper right hand corners of the 
panels. 
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Fig. 2. — The left panel shows the allowed range of 7 for different values of /, regardless of 
the value of To. The right panel shows the allowed range of To, regardless of 7. The ranges 
are obtained using transmitted flux statistics only. 
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are consistently reproduced. 

The constrained parameter space for these three values of / is shown in Figure 1. (We 
mention here that due to finite sampling of the parameter space, the boundaries of the 
allowed region are uncertain by an amount 0.1 along the 7 axis and by 1000K along the T 
axis; this error budget is indicated by a cross at the right top corner of the panels.) It is 
obvious that as we go to lower values of /, the observations allow lower values of T and 
higher values of 7. For example, 1.4 < 7 < 2.3, 0.8 x 10 4 K< T < 1.2 x 10 4 K for / = 0.023 2 
whereas 1.3 < 7 < 1.7, 1.8 x 10 4 K< T < 2.5 x 10 4 K for / = 0.029 2 . It is also seen that 
the area of the allowed region is maximum for / = 0.026 2 and is smaller for higher or lower 
values of /. We mention here that the allowed region is practically zero for / < 0.020 2 and 
/ > 0.032 2 . 

The limits on 7 and T for / in the range 0.023 2 -0.029 2 are shown in Figure 2. The left 
panel shows the allowed range of 7 regardless of the value of T , the right panel shows that 
for T , regardless of 7. 

Although the observations allow / in the range 0.020 2 to 0.032 2 , we find that the match 
between simulations and observations is best for / ~ 0.026 2 . We calculated the \ 2 °f 
the three statistics for different parameter values and found that they are comparatively 
lower for / = 0.026 2 than for higher or lower values of /. Hence, in what follows, we 
shall concentrate on f — 0.026 2 , and see explicitly how all the statistics compare with 
observations. In Figure 3, we show the comparison between simulations and observations for 
various transmitted flux statistics for some particular values of 7 and T . The point to be 
noted here is that our simulations are able to match the observations for all the three statistics 
for a particular range of parameter values. The results obtained using hydrodynamical 
simulations fail to match the observations for the flux correlation function and the power 
spectrum simultaneously (McDonald et al. 2000a), mainly because of the lack of power at 
large scales (due to limited box size). As discussed earlier, our semi analytic simulations 
probe large box sizes without compromising on the resolution, and hence we are able to 
match both the statistics simultaneously. We mention here that the the typical length scales 
probed by both the correlation function and the power spectrum is about 100 h~ x Kpc to 
25 h- 1 Mpc. 

The left panel shows the limit on T for 7 = 1.5. It is clear that for T > 2.2 x 10 4 K, the 
value of F is larger than what is allowed by the observations. At temperatures higher than 
this, the recombination rate is so low that the neutral fraction of hydrogen reduces and hence 
the transmitted flux goes above the allowed limit. Furthermore, we can see that the power 
spectrum also restricts the allowed range of T between (1.2-2.2) x 10 4 K. At higher (lower) 
temperatures, the power at smaller scales are reduced (enhanced) due to excess (less) Voigt 
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f=0.026 2 ,7=l-5 f=0.O26 2 ,T o =1.5xl0 4 K 





Fig. 3. — The comparison between simulations and observed results for / = 0.026 2 . The 
points with errorbars are the observed data points (McDonald et al. 2000a). In the left 
panel, we show the limits on To for a particular 7 (in this case, 1.5), and the right panel 
shows the limit on 7 for Tq = 1.5 x 10 4 K. 
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profile smoothing. ^From the way it is defined (see the previous section), the normalisation 
of the correlation function depends on T . The correlation curves go up when T is decreased. 
However, since the errorbars on £ are comparatively larger, the correlation curve does not 
impose any further constraints. 

The right panel shows the limit on 7 for T = 1.5 x 10 4 K. (We mention here again that 
while changing the value of 7, we change the value of T m also, such that the combination 
7T m , and hence the Jeans length remains unchanged.) The effect of increasing 7 is to 
increase the range of temperature in the IGM for a given baryon density range (Paper I). 
This actually reduces the range in the recombination rate and hence the range of neutral 
hydrogen density. Since a large 7 means less fluctuations in the neutral hydrogen densities, 
there are less number points having extreme values of the flux, as one can see from the PDF. 
It is clear that one can rule out 7 > 2.1 from the PDF. Also, one can see from the flux 
power spectrum curve that there is a reduction (enhancement) in the small scale power for 
larger (smaller) values of 7. This restricts 7 between 1.3-2.1. The correlation curve is quite 
insensitive to 7 as compared to the other two parameters, and hence it does not impose any 
further constraints. Since the normalisation of the correlation function depends on / and To 
but not on 7, the correlation curves are comparatively less sensitive to 7. 

In the allowed ranges of parameters, the match between observations and our simulations 
is quite good for all the three statistics obtained from the transmitted flux. We did not 
compare the simulated flux power spectrum with observations for smaller scales (k > 0.2 
km s" 1 ). The reason for this is the presence of narrow metal lines in the observed spectra, 
which contribute to the small scale power. Detailed discussion regarding this aspect can be 
found in McDonald et al. (2000a). 

Once we have constrained the range of 7 — T space for / = 0.026 2 , it is worth checking 
whether we can match the observed statistics obtained from the Voigt profile decomposition 
of the spectrum. We have used the standard Voigt profile routine (Khare et al. 1997) 
to decompose the observed spectrum into clouds. The minimum number of components 
required to fit an absorption line is constrained by the x 2 minimisation. For this purpose, 
we concentrate on a particular value of T = 1.5 x 10 4 K and 7 = 1.7. For obtaining the 
statistics, we take lines centered around z = 2.26, so as to mimic the observed data of Kim 
et al. (1997) at z = 2.31. The redshift interval considered is Az = 0.26. Figure 4 shows the 
comparison between observations and simulations for the b distribution. The mean b of the 
simulated distribution is 35.19 km s _1 , whereas that of the observed distribution is 36.35 
km s -1 . We have performed a x 2 statistics for the two distributions, and found x 2 jv = 0.61 
(with v = 34, 96.3 per cent likelihood). 

We perform the same exercise with the column density distribution f(Nni)- The com- 
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Fig. 4. — Comparison of the observed b distribution with the simulations. Here, 7 = 1.7, 
T m = 3.01 x 10 4 K, / = 0.026 2 and T = 1.5 x 10 4 K. The observed data is for z = 2.31, taken 
from Kim et al. (1997). 
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parison between simulations and observations is shown in Figure 5. Usually, one assumes a 
power law distribution for f(N HI ), i.e., f(N HI ) oc N H ^ HI . In our case, we obtain the slope 
by carrying out a maximum likelihood analysis (Srianand & Khare 1996) where effects of 
binning are avoided. We present the distribution in the figure in the binned form for the 
purpose of better visualisation. The slope of the column density distribution, (3 hi, in the col- 
umn density range 12.8 < log(A/#//cm~ 2 ) < 16.0 obtained from simulations is 1.31 ± 0.13; 
the corresponding quantity obtained from observations is around 1.35 (Kim et al. 1997). 
One can also compare the values of dN/ dz obtained from simulations and observations. For 
13.77 < log(A^Hi/cm" 2 ) < 16.0, we get dN/dz = 141.06 ± 23.19. The corresponding number 
obtained from observations is between 63.09 and 100.00 (see Figure 2 of Kim et al. 1997). 
For 13.1 < log(iVHi/cm- 2 ) < 14.0, we obtain dN/dz = 202.06 ± 27.75, which is well within 
the observed limits of 158.49 and 223.87. 

The correlation function for the clouds £ c ioud(Av) obtained from our simulations is shown 
in Figure 6, for two different column density thresholds. We have used a velocity bin of 50 
km s" 1 . We have also marked the la and 2a significance levels in the figure, obtained using 
a Poisson distribution. We can see that there is virtually no correlation above 2a significance 
level when the column density threshold is low (A/hi > 12.8 cm^ 2 ). There is a clear positive 
correlation (2.78<r significance) in the velocity bin around Av = 125 km s -1 for clouds with 
A/hi > 10 13 - 8 cm~ 2 . The dependence of clustering on the strength of the lines was noted by 
Cristiani et al. (1997) and Srianand (1997). Kim et al. (1997) find a positive correlation 
of 2.8a significance at the velocity bin 50-100 km s -1 in the observed data. Cristiani et al. 
(1997) too find a positive correlation of about 7a significance at Av = 100 km s _1 , in a wide 
redshift range 1.7 < z < 3.1 using a much larger number of samples. 

While concentrating on the parameter values / = 0.026 2 and T = 1.5 x 10 4 K, we would 
like to see the effect of 7 on the statistics obtained from the Voigt profile decomposition. 
The comparison of the 6-distribution for different values of 7 is given in Table 1. We have 
performed a x 2 test, the results being shown in the same table. The mean b increases with 
7, which is due to the fact that the range of temperatures is higher for large 7. 

Table 1: Comparison between simulated b distribution and observations for / = 0.026 2 , 
T = 1.5 x 10 4 K and different values of 7. The observed mean value of b is 36.35 km s _1 at 



z = 2.31. 
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Fig. 5. — Comparison of the observed column density distribution with the simulations. The 
points in the figure are the observed data. The error bars indicate the values obtained from 
simulations. Here, 7 = 1.7, T m = 3.01 x 10 4 K, / = 0.026 2 and T = 1.5 x 10 4 K. 
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Fig. 6. — The correlation function for clouds obtained from simulations with 7 = 1.7, 
T m = 3.01 x 10 4 K, / = 0.026 2 and T = 1.5 x 10 4 K. The results are presented for two 
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Next, we compare the simulated column density distribution, /(A/hi) with observations 
for different values of 7. Table 2 gives the values of the slope of the distribution, Phi in the 
column density range 12.8 < log(A hi /cm" 2 ) < 16.0 for different 7. The distribution becomes 
steeper as we increase 7. This is consistent with semi analytic results of Paper I. It is very 
clear that the slope for 7 < 1.5 (7 > 2.1) is too flat (steep) to match the observations, even 
within error limits. This is somewhat consistent with what we find from the 6-distribution 
above. We have checked and found that the dN/dz for 13.77 < \og(N HI /cm~ 2 ) < 16.0 for 
different values of 7 is well within the observed range. 

Finally, we discuss the correlation functions. Table 3 shows the correlation function 
within 100 km s" 1 < Av < 150 km s" 1 for different values of 7. There is a slight increase 
in the correlation amplitude as one increases 7. One can, in principle, use this trend to 
constrain the value of 7 through correlation function. However, here we cannot do so be- 
cause of the large errors (a ~ 0.6). We have taken the same number of lines as is done 
in the observations (about 100-140). Consequently, the errors are large and the correlation 
functions are consistent with observations for a wide range of parameter values. 

To summarise, we have shown that, for some particular parameter range, our model is 
consistent with all of the observations (within error limits) obtained from the transmitted 
flux and from the Voigt profile decomposition of the observed spectrum. This justifies our 
approach of modelling the IGM using the lognormal approximation. We have also shown 
that it is possible to put stringent limits on the 7 — T plane for a given / using transmitted 
flux statistics only. 

6. Discussions and Summary 

We have performed a simulation of the Lya absorption spectrum originating from the low 
density IGM using a semi analytic ansatz. We have studied the effect of various parameters 

Table 2: Comparison between simulated column density distribution and observations. The 
observed value of Phi is around 1.35 at z — 2.31 in the column density range 12.8 < 
log(A^ 7 /cm- 2 ) < 16.0 (Kim et al 1997). The value of / and T are 0.026 2 and 1.5 x 10 4 K 
respectively. 
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on the spectrum and the concerned statistics. We have found that our simulations match 
most of the observations available for a narrow parameter range. 

(i) Various statistics performed on the simulated data and the observed points provided 
by McDonald (2000a) over a redshift range 2.09—2.67, constrain the value of / within 0.020 2 - 
0.032 2 , independent of T and 7. In this range of /, we considered three particular values of 
/, namely, 0.023 2 , 0.026 2 and 0.029 2 . We constrain T within (0.8-2.5) xl0 4 K and 7 within 
1.3-2.3. If the range in / is narrowed down through some other studies, the values of 7 and 
To can be constrained further. Although the observations allow / in the range 0.020 2 to 
0.032 2 , we find that the match between simulations and observations is best for / ~ 0.026 2 . 

The values of To and 7 are usually obtained (in previous attempts) from observational 
data through the Voigt profile fitting and the lower envelope TVhi — b scatter plot. The range 
obtained by us is consistent with the one obtained by Schaye et al. (2000). They infer 
1.26 x 10 4 K< T < 2.00 x 10 4 K and 7 = 1.45 - 1.65 for the spectrum of QSO Q1442 at 
z = 2.5 (see their Figure 6). McDonald et al. (2000b) use the lower cutoff of the A H i — b 
scatter plot to infer T and 7. For z = 2.4, they find T = (1.74±0.19) x 10 4 K, 7 = 1.52±0.14 
or T = (1.92 ± 0.2) x 10 4 K,7 = 1.51 ± 0.14, depending on whether they calibrate the data 
using the output from hydrodynamical simulations at z = 3 or z = 2, respectively. On the 
other hand, the hydrodynamical simulations (McDonald et al. 2000a) give slightly lower 
values of To, i.e., To = 1.31 x 10 4 K and 1.6 x 10 4 K, for z = 2 and z — 3, respectively. Ricotti 
et al. (2000) also use the Am — b scatter plot to constrain T between (1-2.4) xl0 4 K at 
z = 1.90 and between (2-2.7) xl0 4 K at z — 2.75. The corresponding constraints on 7 are 
1.32 ± 0.30 at z — 1.90 and 1.22 ± 0.10 at z — 2.75. For clarity, we show how our results 
compare with those obtained by other people in Figure 7. It should be clear from the figure 
and from the above discussion that our results are in quite good agreement with others. 

It should, however, be stressed that the results obtained using the Voigt profile decom- 
position and Nm — b scatter plot have certain inherent biases when compared with those 
obtained from the transmitted flux. For low column density clouds, the error is introduced 
in the b values because of the noise in the spectrum. Also considerable fraction of low col- 
Table 3: The correlation function for Lya clouds within 100 km s _1 < Av < 150 km s _1 for 
/ = 0.026 2 , T = 1.5 x 10 4 K and different values of 7, obtained from the simulations. 
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Fig. 7. — Comparison of the bounds on 7 and T obtained by us with other results. The 
shaded region shows the bound from this work for / = 0.026 2 . The open and the filled squares 
denote the parameter values obtained by McDonald et al. (2000b) depending on whether 
they use the output from hydrodynamical simulations at z = 3 or z = 2, respectively. The 
value obtained by Schaye et al. (2000) for QSO Q1442 around z = 2.5 is denoted by the 
open circle. 
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umn density low b lines that may not trace the low density regions are artificially introduced 
to get a better x 2 while fitting the blends of saturated lines. For the high column density 
clouds, due to saturation, there is a degeneracy between velocity dispersion and number of 
components to be fitted. Our constraints on 7 and T using the transmitted flux statistics 
are free from the above mentioned effects. These constraints on 7 and To can be used si- 
multaneously to constrain the reionisation epoch and the reionisation temperature (Hui & 
Gnedin 1997). 

(ii) We constrain / = (fi^/i 2 ) 2 / J-12 to be in the range 0.020 2 -0.032 2 regardless of the 
values of T and 7. The values for / found by McDonald et al. (2000a) are (0.0257±0.0017) 2 
(for T = 1.31 x 10 4 K) and (0.0239 ± 0.0016) 2 (for T = 1.6 x 10 4 K) for z = 2 and z = 3, 
respectively. This is consistent with the range found in our study. 

The constraint we have obtained on / is important because of the bound it implies on the 
baryon fraction of the universe. The situation is illustrated in Figure 8, where we plot fish 2 as 
a function of J-n- The lower horizontal band corresponds to 0.0170 < VL B h 2 < 0.0208, which 
is considered to be the acceptable range of values from Big Bang Nucleosynthesis (BBN) 
(Buries, Nollett & Turner 2000). As has been noted by several authors and emphasised by 
Padmanabhan & Sethi (2000), this is already in contradiction with the 95 per cent confidence 
limits on Vl B h 2 arising from the analysis of the initial BOOMERANG and MAXIMA data 
The latter bound (0.025 < Q B h 2 < 0.035) is shown in the upper horizontal band in Figure 8 
(for details see Bond et al. 2000; Padmanabhan & Sethi 2000). The bound on Q B h 2 arising 
from the current work (0.020 < Q B h 2 / a/ J-12 < 0.032) is shown as a function of J-12 by the 
curved band running from left bottom to the right top. It is clear from the figure that if 
J_ 12 > 1.2 (indicated by the vertical dashed line in the figure), we have Vt B h 2 > 0.022, which 
is in violation of BBN value. Haardt & Madau (1996), using the QSO luminosity function, 
have estimated J-12 = 1.63 for Q m = 0.2 open universe and J-12 = 1.13 for Q m = 1.0 flat 
universe at z — 2.41. This can be considered as a strict lower bound on J-12, as galaxies 
also contribute equally to the ionising UV background at these redshifts (Steidel, Pettini 
& Adelberger 2000). The bounds obtained from the proximity effect are 0.9 < J-12 < 3.1 
(Scott et al. 2000). Thus, it appears that the bounds on fish 2 obtained from the Lya forest 
analysis could possibly be inconsistent with those obtained from the BBN. 

While this paper was being refereed, three groups have released further data (Netterfield 
et al. 2001,Pryke et al. 2001 and Stompor et al. 2001) with some inital analysis of their 
implications. The BOOMERANG group has given the best bet values of Q B h 2 ~ 0.02 which 
is consistent with BBN results (Netterfield et al. 2001). However, other group still obtains 
Vt B h 2 0.03 (Stompor et al. 2001). There has also been a suggestion (Pettini& Bowen, 
2001) that the a reanalysis of deutrium abundance might raise the BBN bound upwards to 
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Fig. 8. — Comparison of the bounds on fish 2 obtained from our simulations with those 
obtained from Big Bang Nucleosynthesis (BBN) and CMBR analyses. The lower horizontal 
bound shows the region allowed by BBN (Buries, Nollett & Turner 2000) while the upper 
band shows that allowed by initial BOOMERANG and MAXIMA data (Bond et al. 2000; 
Padmanabhan & Sethi 2000). The bound on fish 2 arising from the current work is shown 
as a function of J_i 2 by the curved band running from left bottom to the right top. 



-27- 



about Vt B h 2 ~ 0.025. The situation is therefore unclear at present and these aspects must 
be kept in mind while assessing the importance of the results in Fig. 8. As far as the studies 
on IGM are concerned, we believe that it is important to estimate the value of J_i2 more 
rigorously so as to put a strong constraint on fl B h 2 . Further work in this direction is in 
progress. 
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